function [massProp, flightTime, thrustTime, eclipseTime, coastTime,position] = Low_Thrust_main(mass_sc,T,isp,COE_f,COE_o,eta)

%% Inputs
Tmax    = T;                        % Thrust by thruster
Isp     = isp;                     % Isp of thruster
Mmax    = mass_sc;                     % mass of spacecraft
% thruster_power      = 1350;             %[W] input power in Watt
% thruster_voltage    = 300;              %[V] operation voltage
% thruster_current    = 4.5;              %[A] Operation current in Amps
% battery_capacity    = 3000;         %[W-h] battery capacity
% battery_voltage     = 1500;         %[V] battery provided voltage
UTC_o = [2020 10 01 00 00 00];      % orbital maneuver start time



%% gravitational constant and spacecraft
%--------------- Earth ------------------%
% mu      = 3.986004418*10^5;       % km^3/s^2
% g0      = 9.80665;                % m/s^2
% Re      = 6378.1;   % km
% J2      = 0;
%--------------- Sun ------------------%
mu      = 1.32712440018*10^11;    % km^3/s^2
Re      = 149597870.7;            % km
g0      = 9.80665;                % m/s^2
J2      = 0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

AU2km   = Re;
% Initial orbit (GTO)
a_c   = COE_o(1)/AU2km;                % semi-major axis
ecc_c = COE_o(2);                       % eccentricity
inc_c = COE_o(3)*pi/180;                  % Inclination
ome_c = COE_o(4)*pi/180;                          % Argument of periapsis
Ome_c = COE_o(5)*pi/180;                          % Right ascension ascending node
nu_c  = COE_o(6)*pi/180;                          % True anomaly

% Final orbit (GEO)
a_t   = COE_f(1)/AU2km;
ecc_t = COE_f(2);
inc_t = COE_f(3)*pi/180;
ome_t = COE_f(4)*pi/180;
Ome_t = COE_f(5)*pi/180;
nu_t  = COE_f(6)*pi/180;

TU2sec  = sqrt(AU2km^3 / mu);
TU2days = TU2sec / (60*60*24);
MUkg    = Mmax;

mu   = mu   / (AU2km^3/TU2sec^2);
Re   = Re   / AU2km;
Tmax = Tmax / (MUkg * AU2km * 1000 / TU2sec^2);
Mmax = Mmax / MUkg;
Isp  = Isp  / TU2sec;
g0   = g0   / (AU2km*1000 / TU2sec^2);
c    = Isp * g0;
F    = Tmax / Mmax;

auxdata.mu      = mu;
auxdata.Re      = Re;
auxdata.Tmax    = Tmax;
auxdata.Mmax    = Mmax;
auxdata.Isp     = Isp;
auxdata.g0      = g0;
auxdata.c       = c;
auxdata.F       = F;
auxdata.AU2km   = AU2km;
auxdata.TU2sec  = TU2sec;
auxdata.TU2days = TU2days;
auxdata.MUkg    = MUkg;
auxdata.J2      = J2;

%% Weight definition
Wa = 1;
Wf = 1;
Wg = 1;
Wh = 1;
Wk = 1;

W = [Wa; Wf; Wg; Wh; Wk] / norm([Wa; Wf; Wg; Wh; Wk]);

%% chaser initial state
M0    = Mmax;
[pp_c,ff_c,gg_c,hh_c,kk_c,LL_c] = lowThrustCoe2Mee(a_c,ecc_c,inc_c,Ome_c,ome_c,nu_c);

sc = [pp_c; ff_c; gg_c; hh_c; kk_c; LL_c; M0];


[pp_t,ff_t,gg_t,hh_t,kk_t,LL_t] = lowThrustCoe2Mee(a_t,ecc_t,inc_t,Ome_t,ome_t,nu_t);

st = [pp_t; ff_t; gg_t; hh_t; kk_t; LL_t; M0];

%% Q-Law implementation

del_pp = min((pp_c-pp_t)/10,1/AU2km);
del_ff = 0.0002;
del_gg = 0.0002;
del_hh = 0.0002;
del_kk = 0.0002;
del_cond = [del_pp; del_ff; del_gg; del_hh; del_kk];

sc_history(:,1) = sc;
st_history(:,1) = st;
tt(1) = 0;
t_vec(1) = 0;
Q = proximity_quotient(sc, st, F, mu, W);
Q_history(1) = Q;
cutoff_abs_save = [];
cutoff_rel_save = [];
i = 1;

dL = 10 * pi/180;
tol = 0.1;
tic
count = 0;
while 1

    F = Tmax / sc(7);
    i = i+1;

    %% Q gradient and control definition
    %     [dQdt, u_opt] = proximity_quotient_dot_cutoff(sc, st, F, mu, W);
    [dQdt, u_opt,cutoff_abs,cutoff_rel] = proximity_quotient_dot_cutoff(sc, st, F, mu, W);

    cutoff_abs_save = [cutoff_abs_save,cutoff_abs];
    cutoff_rel_save = [cutoff_rel_save,cutoff_rel];
    if cutoff_rel < eta
        uu = u_opt / norm(u_opt) * 0;
        dQdt = dQdt / norm(u_opt) * 0;
        isthrust = 0;
    else
        uu = u_opt / norm(u_opt) * F;
        dQdt = dQdt / norm(u_opt) * F;
        isthrust = 1;
    end

    %% dt definition
    pp = sc(1);
    ff = sc(2);
    gg = sc(3);
    hh = sc(4);
    kk = sc(5);
    LL = sc(6);

    [~,ecc,~,~,~,nu] = lowThrustMee2Coe(pp,ff,gg,hh,kk,LL);

    rr = pp / ( 1 + ecc * cos(nu) );  % radius from the central body

    dt = dL * sqrt(rr^3/mu);
    %     tt(i) = tt(i-1) + dt;
    %     t_vec(i) = sum(tt);
    if Q + dQdt*dt < 0
        dt = - Q / dQdt;
    end
    tt(i) = tt(i-1) + dt;
    %% chaser trajectory
    k1 = x_dot_efficiency(sc, uu, auxdata,isthrust);
    k2 = x_dot_efficiency(sc + k1*0.5 * dt, uu, auxdata,isthrust);
    k3 = x_dot_efficiency(sc + k2*0.5 * dt, uu, auxdata,isthrust);
    k4 = x_dot_efficiency(sc + k3 * dt, uu, auxdata,isthrust);

    sc = sc + dt*(k1 + 2*(k2 + k3) + k4)/6;
    sc = real(sc);
    sc_history(:,i) = sc;
    %% target trajectory
    k1 = x_dot(st, [0;0;0], auxdata);
    k2 = x_dot(st + k1*0.5 * dt, [0;0;0], auxdata);
    k3 = x_dot(st + k2*0.5 * dt, [0;0;0], auxdata);
    k4 = x_dot(st + k3 * dt, [0;0;0], auxdata);

    st = st + dt*(k1 + 2*(k2 + k3) + k4)/6;
    st = real(st);
    st_history(:,i) = st;
    Q = proximity_quotient(sc, st, F, mu, W);
    Q_history(i) = Q
    if (Q_history(i)-Q_history(i-1))~=0
        count = 0;
    end
    if Q < tol
        break
    elseif (Q_history(i)-Q_history(i-1))==0
        count = count + 1;
        if count == 100
%             fprintf('Q-law can not converge \n\n')
            break
        end
    end
    
end
toc
%% chaser orbit definition
pp_c = ones(1,100) * pp_c;
ff_c = ones(1,100) * ff_c;
gg_c = ones(1,100) * gg_c;
hh_c = ones(1,100) * hh_c;
kk_c = ones(1,100) * kk_c;
LL_c = linspace(LL_c,LL_c+2*pi,100);

[a_c,ecc_c,inc_c,Ome_c,ome_c,nu_c] = lowThrustMee2Coe(pp_c,ff_c,gg_c,hh_c,kk_c,LL_c);
[rr_c,vv_c] = orb2rv(pp_c,ecc_c,inc_c,Ome_c,ome_c,nu_c,LL_c,nu_c+ome_c,Ome_c+ome_c,mu);

%% target orbit definition
pp_t = ones(1,100) * pp_t;
ff_t = ones(1,100) * ff_t;
gg_t = ones(1,100) * gg_t;
hh_t = ones(1,100) * hh_t;
kk_t = ones(1,100) * kk_t;
LL_t = linspace(LL_t,LL_t+2*pi,100);

[a_t,ecc_t,inc_t,Ome_t,ome_t,nu_t] = lowThrustMee2Coe(pp_t,ff_t,gg_t,hh_t,kk_t,LL_t);
[rr_t,vv_t] = orb2rv(pp_t,ecc_t,inc_t,Ome_t,ome_t,nu_t,LL_t,nu_t+ome_t,Ome_t+ome_t,mu);

%% orbit transfer trajectory
pp = sc_history(1,:);
ff = sc_history(2,:);
gg = sc_history(3,:);
hh = sc_history(4,:);
kk = sc_history(5,:);
LL = sc_history(6,:);


[aa,ecc,inc,Ome,ome,nu] = lowThrustMee2Coe(pp, ff, gg, hh, kk ,LL);

[rr,vv] = orb2rv(pp,ecc,inc,Ome,ome,nu,LL,nu+ome,Ome+ome,mu);
rr_q_law = rr;
time_new = linspace(tt(1),tt(end),10000);
time_int = time_new(2);
cutoff_abs_save_new = interp1(tt(1:end-1),cutoff_abs_save',time_new);
cutoff_rel_save_new = interp1(tt(1:end-1),cutoff_rel_save',time_new);

rr_new = interp1(tt,rr',time_new);
rr = rr_new';

%% Find eclipse information
semiMajor = 6378.137;
eclipse_save = [];
for ii = 1:length(time_new)
    MJD_TT = JulianDate([UTC_o(1:5),UTC_o(6)+time_new(ii)]);
    r_sun = Sun(MJD_TT);
    r_save(:,ii) = r_sun;
    R_eci = rr(:,ii)*semiMajor;
    eclipseCondition1 = dot(R_eci,r_sun)/...
        (norm(R_eci)*norm(r_sun));
    eclipseCondition2 = norm(R_eci)*sqrt(1-eclipseCondition1^2);
    if (eclipseCondition1<0)&&(eclipseCondition2<semiMajor)
        eclipse_save = [eclipse_save,ii];
    end
end

%% trajectory plot
coast_arc_ind = find(cutoff_abs_save_new < eta);
thrust_arc_ind = find(cutoff_abs_save_new >= eta);

%%
position.coast=rr*AU2km;
position.thrust = rr(1:3,thrust_arc_ind)*AU2km;
massProp = MUkg-sc_history(7,end)*MUkg;
flightTime = tt(end)*TU2days;
thrustTime = time_int*length(thrust_arc_ind)*TU2days;
coastTime = time_int*length(coast_arc_ind)*TU2days;
eclipseTime = time_int*length(eclipse_save)*TU2days;
% figure;
% plot(1:length(Q_history),Q_history)
%
propellant_mass = MUkg-sc_history(7,end)*MUkg;
% fprintf('Necessary Propellant mass: %5.2f Kg\n\n',propellant_mass);
% fprintf('Total Flight Time: %3.2f Days\n\n',tt(end)*TU2days);
% fprintf('Total Thrusting Time: %3.2f Days\n\n',time_int*length(thrust_arc_ind)*TU2days);
% fprintf('Total Coasting Time: %3.2f Days\n\n',time_int*length(coast_arc_ind)*TU2days);
% fprintf('Total Eclipse Time: %3.2f Days\n\n',time_int*length(eclipse_save)*TU2days);